This exercise simulates the motion of rigid bodies, and their interactions through collision. The objectives are:
The repository provides a skeleton code on which to build your first exercise. CMake allows you to work and submit your code in all platforms. The entire environment is in C++, for most part, you only code the mathematical-physical parts.
The basic scenario is limited to convex objects, and includes no air resistance or friction. The environment loads scenes that describe a set of objects, and their positions in space, and the world contains a big, immobile, and heavy plate on which the objects fall. There are no ambient forces in the basic setting other than gravity.
Collisions are limited to a single point per object. (If more, one point is chosen -- this might lead to mildly non-physical behavior). The environment is already configured to run in a time loop, where it integrates movement and detects collisions in each step.
See below for details on where to do all that in the code.
The above will earn you 70% of the grade. To get a full 100%, you must choose and augment the implementation by 2 of 6 extension options below. Some will require minor adaptations to the GUI or the function structure. Each extension earns up to 15%. The exact grading will be commensurate with the difficulty. Note that this means that all extensions are equal in grade; if you take on a hard extension, it's your own challenge to complete well.
You may invent your own extension as substitute to one in the list above, but it needs approval by the Lecturer.
The skeleton uses the following dependencies: libigl, and consequently Eigen, for the representation and viewing of geometry, and libccd for collision detection. libigl viewer uses nanogui for the menu. Everything is bundled as either submodules, or just incorporated code within the environment, and you do not have to take care of any installation details.
To download the library, use:
git clone --recursive https://github.com/avaxman/INFOMGP-Practical1.git
To compile the environment, go into the practical1
folder and enter in a terminal (macOS/Linux):
mkdir build
cd build
cmake -DCMAKE_BUILD_TYPE=Release ../
make
In windows, you need to use cmake-gui. Pressing twice configure
and then generate
will generate a Visual Studio solution in which you can work. The active soution should be practical1_bin
. Note: it only seems to work in 64-bit mode. 32-bit mode might give alignment errors.
You do not need to utilize any dependency on your own, or install anything other than the above. For the most part, the dependencies are parts of code that are background, or collision detection code, which is not a direct part of the practical. The exception is Eigen for the representation and manipulation of vectors and matrices, but it is a quite a shallow learning curve. It is even possible to learn most necessary aspects from looking at the existing code. However, it is advised to go through the "getting started" section on the Eigen website (reading up to and including "Dense matrix and array manipulation" should be enough).
All the code you need to update is in the practical1
folder. Please do not attempt to commit any changes to the repository. You may ONLY fork the repository for your convenience and work on it if you can somehow make the forked repository PRIVATE afterwards. Solutions published online disqualify!
Most of the action happens in scene.h
. The main function is:
void updateScene(double timeStep, double CRCoeff, MatrixXd& fullV, MatrixXi& fullT){
fullV.conservativeResize(numFullV,3);
fullT.conservativeResize(numFullT,3);
int currVIndex=0, currFIndex=0;
//integrating velocity, position and orientation from forces and previous states
for (int i=0;i<rigidObjects.size();i++)
rigidObjects[i].integrate(timeStep);
//detecting and handling collisions when found
//This is done exhaustively: checking every two objects in the scene.
double depth;
RowVector3d contactNormal, penPosition;
for (int i=0;i<rigidObjects.size();i++)
for (int j=i+1;j<rigidObjects.size();j++)
if (rigidObjects[i].isCollide(rigidObjects[j],depth, contactNormal, penPosition))
handleCollision(rigidObjects[i], rigidObjects[j],depth, contactNormal.normalized(), penPosition,CRCoeff);
//Code for updating visualization meshes
for (int i=0;i<rigidObjects.size();i++){
fullT.block(currFIndex, 0, rigidObjects[i].T.rows(),3)=rigidObjects[i].T.array()+currVIndex; //need to advance the indices, because every object is indexed independently
fullV.block(currVIndex, 0, rigidObjects[i].currV.rows(),3)=rigidObjects[i].currV;
currFIndex+=rigidObjects[i].T.rows();
currVIndex+=rigidObjects[i].currV.rows();
}
currTime+=timeStep;
}
The two most important functions are integrate()
and handleCollision()
. They all contain a mixture of written code, and code you have to complete. The code you have to complete is always marked as:
/***************
TODO
***************/
For the most part, the description of the function will tell you what exactly you need to put in.
The program is loaded by giving a TXT file that describes the scene as an argument to the executable. The file should be in the data
subfolder, which is automatically discovered by the CMake. The format of the file is:
#num_objects
object1.off density1 is_fixed1 COM1 o1
object2.off density2 is_fixed2 COM2 o2
.....
Where:
objectX.off
- an OFF file (automatically assumed in the data
subfolder) describing the geometry of the a triangulated mesh. Meshlab can view these files, but their format is pretty straightforward. The object is loaded, and translated to have a center of mass of (0, 0, 0) in its object coordinates. density
- the uniform density of the object. The program will automatically compute the mass by the volume. is_fixed
- if the object should be immobile (fixed in space) or not. COM
- the initial position in space where the object would be translated to. That means, where the COM is at time t = 0. o
- the initial orientation of the object, expressed as a quaternion that rotates the geometry to o * object * inv(o) at time t = 0.The viewer presents the loaded scene, and you may interact with the viewing with the mouse: rotate with the left button pressed and moving around (the "[" and "]" buttons change the behaviour of the trackball), zoom with the mousewheel, and translate with the right button pressed and dragging. Some other options are printed to the output when the program starts.
The menu also controls the visual features, and the setting of the coefficient of restitution and the time step. They can be updated at any point in the simulation. You might add more parameters with some extensions. Everything is set up in main()
.
The simluation can be run in two modes: continuously, toggled with the space
key (to stop/run), and step by step, with the S
key. This behavior is already encoded. The visual update of the scene from the objects is also already encoded in updateScene()
There are two main classes, scene
, and RigidObject
. Both are in scene.h
, and will be updated by you. They are commented throughout, so their individual properties are understood. Each mesh, and the platform, are rigid bodies of their own. The geometry is encoded as follows:
MatrixXd origV; //original vertex positions, where COM=(0.0,0.0,0.0) - never change this!
MatrixXd currV; //current vertex position
origV
and currV
are ∣V∣ × 3 matrices encoding all the vertices of the mesh, row by row, and T encodes the indices to the triangles (you do not need to use it directly). You should never update origV
. currV
should be updated to reflect the result of every time step, and this is what you see on screen.
Quaternions represent orientations and rotations, where if the neutral (initial) orientation of a vector is v, and the orientation quaternion is q, then the final orienation is qvq − 1. The QRot
function in the auxfunctions.h
file implements that (and several other functions for quaternions are available in that file). Note the function Q2RotMatrix
that produces the rotation matrix for that orientation, and should be used for the transformation of the inertia tensor.
You do not have to compute the entire algorithmic environment from scratch. The things that you are given are:
getCOMandInvIT
that computes the original COM and the inverse inertia tensor from the original OFF file, and is called by the RigidObject
constructor. You do not need the COM it computes; the constructor translates the object (origV
coordinates) to the origin, so it always has COM = (0, 0, 0). The constructor also create currV
as a translation and rotation of origV
to fit the prescribed values from the scene file.The inverse inertia tensor you get from getCOMandInvIT
is not after applying the orientation, not even that in the scene file. That is, what you get is the inverse inertia tensor of origV
around the origin. You will have to compute the inverse inertia tensor for a given currV
, according to the its current orientation, and it is always then around the COM of the moving object. Be careful to apply the correct rotation!
Here are detailed answers to common questions. Please read through whenever ou have a problem, since in most cases someone else would have had it as well.
Q: I am getting "alignment" errors when compiling in Windows. A: Delete everything, and re-install using 64-bit configuration in cmake-gui
from a fresh copy. If you find it doesn't work from the box, contact the Lecturer. Do not install other non-related things, or try to alter the cmake.
Q: I have an angular velocity ω⃗, how do I integrate it? A: There are actually two ways to do that that could be considere "forward Euler". The following is the simplest (radial integration), and working in the demo:
QExp
given to you in the code. However, you need to feed it with the imaginary quaternion Δq = QExp((0, Δt ⋅ ω⃗)).QMult
for quaternion multiplication.